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Abstract 

We address the dynamics of adsorbed molecules (a fundamental issue in surface physics) within 
the framework of a Master Equation scheme, and study the diffusion of particles in a finite cubic 
lattice whose boundaries are at the z = 1 and the z = L planes where L = 2,3,4,..., while the 
X and y directions are unbounded. As we are interested in the effective diffusion process at the 
interface z = 1, we calculate analytically the conditional probability for finding the system on the 
z = 1 plane as well as the surface dispersion as a function of time and compare these results with 
Monte Carlo simulations finding an excellent agreement. 

PACS numbers: 05.40.Fb, 02.50.Ey, 05.10.Ln, 46.65.+g 
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I. INTRODUCTION 



The dynamics of 



The dynamics of adsorbed molecules is a fundamental issue in interface science 
and is also crucial in a large number of emerging technologies 



Adsorption at the solid-liquid interfaces arises, for instance, in many biological contexts 
involving protein deposition ^ 13. in solutions or melts of synthetic macromolecules 
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loidal dispersions 19|, and in the manufacture of self-assembly mono- 



and multi-layers (ll|, ^ 

ExpeHmenta. stud.es of surfactant molecules BhO. proteins BMM and synthefc 

polymers p, HO] confined to interfaces have identified different types of surface translational 
motions. One is in-surface self-diffusion of individual molecules, which has been investigated 
with fluorescence recovery after photobleaching (FRAP) methods ji, 10, 2^. Measured self- 
diffusivity at liquid-solid interfaces are much smaller than bulk values [2J] and similar to 
bulk values in liquid-gas cases !l^. A second source of motion, exclusive to the liquid-fluid 
interface, is surface visco-elasticity. Adsorbed surface phases possess compressibility and 




axation kinetics 
, with 



|24, 

25[ providing 



viscosities which govern the dynamics of the surface density waves. Their re 
have been intensively studied experimentally and theoretically Jl| 

direct viscoelastic measurement j5| and surface light scattering studies 
considerable support for current theories. 

In this paper we explore another mechanism called hulk-mediated surface diffusion. This 
mechanism arises at interfaces separating a liquid bulk phase and a second phase which may 
be either solid, liquid, or gaseous. Whenever the adsorbed species are soluble in the liquid 
bulk, adsorption-desorption processes occur continuously. These processes generate surface 
displacement because desorbed molecules undergo Fickian diffusion in the liquid's bulk, and 
are later re-adsorbed elsewhere. When this process is repeated many times, an effective 
diffusion results for the molecules on the surface. The importance of bulk-surface exchange 
in relaxin g ho mogeneous surface density perturbations is experimentally well established 

Here, unlike our previous work [27] on the adsorption-desorption mechanism, we assume 
that the bulk is finite. The main goal of this work is the observation of an effective diffusion 
process at the interface z = 1. We calculate analytically its variance, {r'^{t)), as a function 
of time and P{x, y, z = 1; t\0, 0, 1, t = 0), the conditional probability of being on the surface 
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at time t since the particle arrived at t — 0, that we indicate as P^^iit). 

The organization of the paper is as follows. In the next section we describe the model by 
means of a set of Master Equations and solve this system in Fouricr-Laplacc space. After 
that, we present a special case (the bilayer case) where we can obtain explicit analytical re- 
sults for P^=i(t) and {r^{t)). Next we compare these theoretical results with numeric Monte 
Carlo simulations. Moreover, we present results obtained by numerical inverse transforma- 
tion of the equations for (r^) and Pz=i for the cases with more layers. In the last section we 
discuss the results and draw some conclusions. 

II. THE MODEL 

Let us start with the problem of a particle making a random walk in a finite cubic lattice. 
The bulk is bounded in the z direction where the particles can move from z — 1 to z — L. 
The X and y directions remain unbounded. The position of the walker is defined by a 
vector r whose components are denoted by the integer numbers n, m, I corresponding to the 
directions x, y and z respectively. 

The probability that the walker is at (n, m, I) for time t given it was at (0, 0, Iq) at t = 0, 
P{n, m, I; t\0, 0, Iq, t = 0) = P{n, m, I; t), satisfies the following master equation 

P{n,m^l]t) = jP^n^mJ + l;t) — SP{n,m,l;t) 

+a^[P{n - 1, m, I; t) + P{n + 1, m, I; t) - 2P{n, m, I; t)] 

+/3^[P(n, m - 1, Z; t) + P(n, m + l,l;t)- 2P(n, m, t)], for Z = 1 

P(n, m, /; t) = a[P{n - 1, m, /; t) + P{n + 1, m, /; t) - 2P{n, m, I; t)] 
+l3[P{n, m-l,l;t) + P{n, m + l,l;t)- 2P(ra, m, I; t)] 
+-fP{n, m,l + l;t) + 5P{n, m,l-l;t)- 2-fP{n, m, t), for / = 2 

P(n, m, I; t) = a[P{n - 1, m, /; t) + P{n + 1, m, /; t) - 2P(n, m, /; t)] 
+/3[P(n, m-l,l]t) + P{n, m + l,l;t)- 2P(n, m, I; t)] 

+7[P(n, m,l + l;t) + P{n, m,l-l;t)- 27[P(n, m,l;t)], for 3 < / < L - 1 
P{n,m,l;t) = 7P(n, m, / — 1; t) — 7P(n, m, /; t) 

+a[P{n - 1, m, I; t) + P{n + 1, m, I; t) - 2P(n, m, I; t)] 

+p[P{n,m-l,l;t) + P{n,m + l,l;t) -2P{n,m,l;t)], for / = L. (1) 
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where a, p and 7 are the bulk transition probabihties per unit time in the x,y and z directions 
respectively, and 6 is the desorption probability per unit time from the boundary plane z = 1. 

It is important to note that the model presented in Eq. (P), allows the possibility that the 
particles can move in the plane z = 1 with temporal frequencies in the x direction and (3^ 
in the y direction. If these temporal frequencies are equal to zero, the motion through the 
z = 1 plane is exclusively due to the dynamics across the bulk. In addition we can observe 
that this is a finite set of L equations. This fact establishes an important difference with the 
infinite set of equations presented in [^^l , a crucial difference because this generates different 
solutions to the problem. In order to solve this finite set of equations, we take the Fourier 
transform with respect to the x and y variables, and the Laplace transform in the t variable. 
After these transformations, we obtain the following set of equations 

sG{kx-, ky, I] s) P{kxi ky, /, t = 0) = 'jGi^kx, ky, / + 1; s) 

-6G{kx, ky, I] s) + A^{kx, ky)G{kx, ky, I; s), for / = 1 
sG{kx, ky, I', s) P{kxi ky, l,t = 0) = A(^kx, ky)G(^kx, ky, t, s) + 5G{kx, ky, I 1; s) 

+-iG{kx, ky, l + l;s) - 2-fG{kx, ky, I; s), for / = 2 
sG{kxi ky, I', s) P{kx, ky, l,t = 0) = A(^kx, ky)G{kx, ky, I; s) + '^\G{kx, ky, I 1; s) 

+'^G{kx, ky, l + l;s)- 2-fG{kx, ky, l;s)], for 3 < / < L - 1 
sG{kxi ky, l] s) P{kxi ky, l,t = 0) = A(^kx, ky)G(^kx, ky, I', s) 




for / = L 



(2) 



where we have defined 



G{kx, ky, l] s 



G{kx,ky,l;s\0,0,lo;t = 0) 




00 




(3) 



L indicates the Laplace transform of the quantity inside the brackets, and 




(4) 




(5) 
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Using the matrix formalism, Eq. Q can be written as 

[si - H]G = 5u, = III,, 
where G is an L x L array that has the following components 

Gil, = [G[k^, ky, I] s\n, m, /q; to]], 
/ is the identity matrix and H is a. tri-diagonal matrix whose components are 



(6) 



(7) 
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The G parameter is defined as 
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(8) 



G = -2^ + A{k,,ky). 

In order to find the solution to the Eq. we decompose the H matrix according to 

H = A{k,, ky)i +Hq + H, + H^, (9) 



where 



^^0 



/-27 7 O--- 

7 -27 7 O--- 

^ -27 7 ■ ■ ■ 









\ 











V 



7 —27 7 
7 -7 / 



Ai 



1 ifi = j = 1 

otherwise 



2 



1 if i = 1 and j = 2 
otherwise 



and 



We also define 



Ai = -S - [-2-f + A{k,,ky) - A\K,ky)], 

A2 = 5-7. (10) 

= [si-iA{k^,ky)i+Ho)]-\ 

= [si-iA{k,,ky)i+Ho + Hi)]-\ (11) 
A formal solution of the Eq. (jH)) is 

G=[si-H]-\ (12) 

We can show, by applying the Dyson formula |31], that 

A2 Gj^ G\i^ 



1 - A2 G\^ ' 



Gilo - (^lio + — / ^1 ° ' (13) 

'-il _ qO 

nio - iio ■ QO 



GL = Gl + (14) 
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The solution for Gii^ can be obtained analytically |2^. The result is 



L-l 



Gl E /'^/'o.2^ ^ _ ^(^^^ _ 27cos(g.) ' ^^^^ 

where 

fii = Ksm{lqi), (16) 

and 



(2z + 1) 



TT 



» 2LH-1- (1^' 

The constant K is obtained by exploiting the orthonormality relations for the fu functions 
given by 
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while the completeness relation for the fu functions is written as 

J2f^^f^^ = ^^r (19) 

i=l 

The corresponding expression for K is 

We are now able to find the statistical quantities which describe the diffusion problem over 
the surface. We are interested in the probability of finding a particle in the site (m, n, / = 1) 
at time t given it was in (0,0,/ = 1) at t = 0. This quantity is obtained by applying the 
inverse Laplace transform on the Gu matrix element. Another direct measurable experi- 
mental magnitude is the variance {{r^{t))) of the probability distribution at time t 
over the plane 

{r\t))piane. (21) 

This quantity measures the particles dispersion over the surface and has a direct relation to 
the diffusion coefficient. When P(m, n, / = l;t|0,0,/o = l;t = 0) is known, the variance is 
calculated in the following manner 

oo 

{r\t))piane= P{m, Ti, I = 1; t\0 , , = t = 0){m^ + u^) . (22) 



m,n=—oo 



In order to obtain this expression, we have assumed symmetric properties for the diffusion 
along the x and y directions, that is < x{t) >=< y{t) >= 0. 
Finally the variance in the Laplace space can be obtained as 

< r^{s) >plane= + Q^Wn] \k:.=ky=0 ■ (23) 

X y 

III. THE BILAYER CASE 

So far, we have developed a general theory describing the important problem of the 
diffusion of a particle system over a surface which is surrounded by a bounded bulk. The 
expressions we found, are in the Laplace space and the obtention of the inverse transform 
is usually a non trivial task. 

In this section we restrict the problem to finding a particular solution, considering the 
case of a system of particles moving inside a bilayer, that is the space formed by the surface 



where we investigate the movement of the particles and a second layer. To obtain the 
P{m,nJ = l;t|0,0,/o = 1;^ = 0) {Pz=i{t)), we take the expression in Eq. (fT^ with L = 2 
and evaluate for = ky = 0. Then we obtain the inverse Laplace transform. The result is 

This expression establishes that this probability is only a function of the temporal adsorption 
and desorption rates and does not depend on the temporal rates on the plane. From this 
equation, we can also obtain the initial condition of the problem by evaluating Eq. (j24|l at 
t = 0, obtaining Pz=i{t = 0) = 1. 

If we consider the long time limit, the expression has the following behavior 

lim P,=i(t) = (25) 

This expression shows that, at this limit, the probability of finding the particle in the z = 1 
plane does not decay to zero (as happens in unbounded systems |23|), but it approaches an 
asymptotic value which is a function of the ratio of the temporal rates. 

If the adsorption rate 7 is large (compared with the desorption rate 5), then the probabil- 
ity is large, the particles go back to the plane frequently. On the other hand, if the desorption 
rate is large, the particles go away from the surface and so the probability decays. These 
behaviors are shown and quantified by Eq. (j23}. 

In the case of the variance, the result we have obtained can be decomposed as follows 



{r\t)) = {r\t))^iane + {r\t)),oU (26) 



where 



275 + ^ ^ . 72 



(7 + 6)- 



! 2 \ 

'-t exp[-(7 + <5)t] + ^^^^tJ, (27) 



{r\t)u = («+/?) {j^^, (^^p [-(7 + m - 1) 

' texp[-i^ + 6)t] + -^^t]. (28) 



From Eqs. (j^Bj) . and ()28|) we can make the following observations. Firstly, the variance 
can be decomposed into two contributions, one corresponding to the particle movement 
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across the bulk and the other to the surface movement. The dependence of the variance on 
the time rates parallel to the surface is linear, while in these same relations, the adsorption 
and desorption rates enter in a more complicated way. The functional form obtained is 
similar for both movements (which however are nonindependent). 

For a large evolution time, the mean square distance or dispersion grows linear with 
time; each diffusive process has its own slope or growing rate (these slopes are associated 
with the diffusion coefficient). However, this is an expected behavior due to the model 
we are using. The other contributions are transient ones, that decay with a time constant 
T = (7 + 5)^^, hence, the stronger the adsorption or desorption, the faster is the decay of 
these contributions. 

IV. RESULTS 

In this section we show the theoretical results, including some special numerical proce- 
dures, and make the comparison with Monte Carlo simulations. In all simulations we have 
fixed the following parameters: a = /3 = 7 = 1 and = f3^ = 0, and have averaged over 
10^ realizations. 

In Fig. H we present the theoretical-numerical (that is using a computer program to 
calculate the inverse Laplace transform j^^l ) , theoretical results obtained from Eq. and 
simulation results for the temporal evolution of the probability to find the system on the 
surface for the bilayer case. Here we show the curves for two values of the desorption rate 
6. As is apparent from the figure, there is excellent agreement between the theoretical and 
simulation results. Such excellent agreement indicates that the numerical procedure for the 
obtention of the inverse Laplace transform is a reliable tool, and that we can trust their 
results in those cases where analytical results are not accessible (for instance the cases with 
larger number of layers that we will consider in the following). In Fig. |2l and again for 
the bilayer case, we depict for the variance ((r^(t))), both theoretical and numerical results. 
Again the agreement is excellent. 

In Figs. 01 and m we present the P{m,nJ = l;t|0,0,/o = 1;^ = 0) and the (r^(t)) 
but now for a number of layers larger than two. Here we compare the theoretical and 
numerical results again. We remark that the theoretical results were obtained fitting the 
inverse Laplace transformation of Eqs. (fTH|) and (??) numerically. Again we have found an 
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FIG. 1: Temporal evolution of the Pz=i{t). We have shown two cases: (i) A represent theoretical 
points (see Eq. (|24() ). the continuous line indicates the theoretical- numerical results and Q ^i'^ the 
simulations data for 6 = 0.1; (ii) T corresponds to theoretical points, the dashed line represents 
theoretical-numerical results and □ the simulations data for 5 = 0.5. 

excellent agreement between theoretical and simulation results. 

Figures El IHl and U\ correspond to an analysis for the asymptotic behavior of the finite 
system. Figure El depicts the curve obtained for the (r^) as a function of k (the number of 
layers) for three different and large observational times: t = 1200, 1300, 1500. The insert 
shows that for these times, the system is well inside the asymptotic region. As can be seen 
from the figure, there is a maximum in the motion of the system. In other words, there is 
an "optimal" number of layers for which the system can spread more rapidly. For a larger 
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FIG. 2: Time evolution of (r^). We have represented two cases: (i) A represent theoretical points 
(see Ea. ()24() ). the continuous line indicates the theoretical- numerical results and Q ^i'^ the simula- 
tion data for 5 = 0.1; (ii) T correspond to theoretical points, the dashed line represents theoretical- 
numerical results and □ the simulations data for 5 = 0.5. 

number of layers, the system converges to an asymptotic limit. This is reasonable due to 
the fact that the finiteness of the system tends to disappear. 

Figure ini depicts the same behavior as Fig. |31but now we have fixed an observational time 
(t = 1500) and we have used the desorption rate (6) as parameter. The figure also shows the 
maximum on the number of layers again but from this figure we can see that this maximum 
moves towards the lower layers as the desorption rate increases. The insert shows the time 
evolution of (r^) for large values of k. We compare this behavior with the one corresponding 
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FIG. 3: Time evolution of the Pz=i{t). We have represented two cases: (i) the continuous hne 
depicts the theoretical result obtained numerically and Q the simulation points for 3 layers; 
(ii) the dashed line represents theoretical results and □ the simulations for 4 layers. 

to the infinite case From this figure we can say that the finite system behaves, for the 
used parameters, as an infinite one when the number of layers is k > 50. 

If the parameter 6 is increased, the maximum finally disappears as can be seen in Fig. [7| 
This figure shows two curves for two different 6 for t = 1500. The insert shows the Pz=i{t) 
for these desorption rates. 

Finally we show a fitting of {r'^{t)) as a function of k for a large evolution time. For 
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FIG. 4: Time evolution of the variance (r^). We have represented two cases: (i) the continuous hne 
depicts the theoretical result obtained numerically and Q the simulation points for 3 layers; 
(ii) the dashed line represents theoretical results and □ the simulation for 4 layers. 

fitting purposes we have used the following function 

(r'it)) = C f, (29) 

where C is a constant and e is the fitting parameter. In Fig. IHlwe depict the results obtained 
for e function of k. 

We observe two regions: for k < 20 the particles describe a diffusion movement. In 
particular this behavior was shown for the bilayer model. The insert of the figure shows 
a zoom for a small number of layers, and for two different times. The linear behavior is 
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FIG. 5: (r^) vs k for the case 6 = 0.02. The Q correspond to an observation time t = 1200, ■ is for 
t = 1300 and + for t = 1500. The insert shows the Pz=i{t) vs k. The system is in the asymptotic 
region. 

apparent. A second region appears for k > 20 where the movement results subdiffusive. For 
large k the e parameter approaches an asymptotic value equal to 0.5. This situation was 
predicted [^l for infinite systems. 

V. CONCLUSIONS 

We have studied the evolution of particles diffusing in a volume, but have analyzed the 
statistical properties of their evolution on a surface. The diffusion can be performed on 
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FIG. 6: (r^) vs k. We have represented the case for t = 1500. The Q correspond to a desorption 
rate 5 = 0.02, while ■ is for 5 = 0.04 and + for 5 = 0.06. Insert: time evolution of (r^) for 5 = 0.02 
and K = 100 (O)) ^-i^d for 5 = 0.06 and k = 100 (□). The continuous lines correspond to the 
infinite bulk behavior [2^ . 

the surface itself or across the bulk surrounding the surface. In this work we focus on the 
particle diffusion due to the movement across the bulk, that is we have considered the hulk 
mediated surface diffusion of the particles. The main feature of the bulk is its finiteness in 
one direction (the axial one). The other directions are unbounded. This work complements 
a previous study [2^ in which we have analyzed the particle diffusion in a semi-infinite bulk. 

Here we have presented a theoretical model based on a set of Master Equations which 
describe the movement of particles over a simple cubic lattice. This model is general in 
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FIG. 7: (r^) vs k, The O corresponds to a, 6 = 1.5, ■ to a 5 = 2. Insert: Pz=i{t) vs k for these 
desorption rates. 

the sense that we include both kind of particle movement: on the surface and on the bulk. 
We solved this problem by using techniques of the Laplace and Fourier transformation and 
exploiting Dyson's formula. It is worth remarking here that the model can describe the 
evolution of the system everjrwhere in the bulk, and the evolution of the system for all time. 
We particularize to the study over a single plane, i.e. the surface. We have obtained general 
solutions in the Laplace space for the Pz=i{t) and the {r'^{t)) for any kind of bulk. 

In order to find an analytical solution in space and time, we have particularized the 
problem to the bilayer case. The expressions obtained with this assumption were compared 
with Monte Carlo simulations, finding excellent agreement. Moreover we were able to handle 
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FIG. 8: e vs k. Circles correspond to 6 = 0.5 and squares to 6 = 2. The time evolution was 
t = 2000. Insert: e vs k for a small number of layers and t = 2000 (■), t = 3000 (Q)- The line is 
a linear fitting of the data. 

cases with more layers. In these cases, we have obtained the inverse Laplace transform 
using numerical methods, and compared these results with Monte Carlo simulations. The 
agreement between both results are excellent again. As indicated above, the numerical 
procedure to obtain the inverse Laplace transform is a reliable tool, and we have shown that 
we can trust their results in those cases where analytical results are not accessible. 

Finally we have studied the behavior of the spreading in the asymptotic region as a 
function of the number of layers of the system. We observed an "optimal" number of 
layers for which {r^{t)) reaches its maximum value. It is worth remarking here that this 
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effect occurs in the limit of "strong adsorption", that is when the relation ^ is small, and 

disappears in the limit of "weak adsorption", that is when ^ is large. ^ ^ 

A possible generalization of the present and the previous related work [22i] consists in 
considering the possibility of non Markovian dynamics, and study the effect of such dynamics 
on the statistical features of the system. This is the subject of further work. 
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